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Abstract — We present a hierarchical computation approach 
for solving finite-time optimal control problems using operator 
splitting methods. The first split is performed over the time 
index and leads to as many subproblems as the length of the 
prediction horizon. Each subproblem is solved in parallel and 
further split into three by separating the objective from the 
equality and inequality constraints respectively, such that an 
analytic solution can be achieved for each subproblem. The 
proposed solution approach leads to a nested decomposition 
scheme, which is highly parallelizable. We present a numerical 
comparison with standard state-of-the-art solvers, and provide 
analytic solutions to several elements of the algorithm, which 
enhances its applicability in fast large-scale applications. 

I. INTRODUCTION 

Online optimization and optimal control methods are in- 
creasingly being considered for fast embedded applications, 
where efficient, reliable, and predictable computations in- 
volved in calculating the optimal solutions are a necessity. 
The potential use of optimal control in such embedded 
systems promises energy savings and more efficient resource 
usage, increased safety, and improved fault detection. The 
range of application areas that can benefit from embedded 
optimization include the mechatronics, automotive, process 
control and aerospace sectors [1]. The promise of unprece- 
dented performance and capabilities in these applications, 
which typically rely on large-volume, real-time embedded 
control systems, has fueled recent research efforts towards 
fast and parallel optimization solvers. 

One of the main research directions aim at develop- 
ing special-purpose optimization solvers that target typical 
control or estimation problems arising in optimal control. 
Parallel solutions to systems of linear equations appearing 
in interior-point, and active set methods have been studied 
in [2]-[5]. In this work we consider a quadratic finite- 
time optimal control problem for discrete-time systems with 
constrained linear dynamics, which appears in typical model 
predictive control problems [6]. We investigate and develop 
different parallelizable algorithms using operator splitting 
techniques [7], [8] that have recently shown great promise for 
speeding up calculations involved in computing optimal solu- 
tions with medium accuracy [9]-[l 1]. Our approach relies on 
a hierarchical splitting up of the specially structured finite- 
time optimal control problem. The first split is performed 
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over the time index and leads to as many subproblems as 
the length of the prediction horizon. Each subproblem can 
then be solved in parallel and further split into three by 
separating the objective from the equality and inequality 
constraints respectively, such that an analytic solution can 
be achieved for each subproblem. The proposed solution 
approach leads to a nested decomposition scheme, which 
is highly parallelizable. The proposed three-set splitting 
method does not only solve the particular quadratic programs 
(QPs) that appear in the update steps of the time- splitting 
algorithm efficiently, but also provides a compact, standalone 
alternative for solving generic QPs. 

The paper is structured as follows. Section UTJ presents 
the main idea behind the time- splitting optimal control 
approach for parallel computations, using the Alternating 
Direction Method of Multipliers and deriving the exact 
formulas required for each subproblem and update step. In 
Section |llll we propose an alternative scheme for solving the 
QPs with general polyhedral constraints that arise in the time- 
splitting update steps (or for any other generic QP). The 
two splitting schemes are combined in a hierarchical fashion 
in Section [IVl and numerical experiments are performed in 
Section [V] to compare its performance with some advanced 
solvers in the literature. Section |VT] concludes the paper. 

II. TIME-SPLITTING OPTIMAL CONTROL 

A. Problem Formulation 

We consider the following finite-time optimal control 
problem formulation that arises in typical model predictive 
control applications: 



1 N / \ 

minimize - ^ [xj Qx t + uj Ru t J 

2 1=0 



(la) 



subject to (x t ,u t )e3£ t x%, t = 0,...,N (lb) 
x t+ \ = A t x t + B t u t + c t ,t = 0, ... ,N - l(lc) 

where the decision variables are the states x t G R", and the 
inputs u t G R m of the system for t =0,...,N. The index t 
denotes time, and the system evolves according to linear 
dynamics constraint (fTcT> where c t G R" is considered to be 
a known disturbance. Here, N is the prediction horizon and 
Q G R nx '\R G R mxm are symmetric matrices. The stage cost 
functions in (flat are convex quadratic with Q^O and R y 0. 
The stage-wise state-input pairs are constrained to reside 
within polyhedra (flbl denoted by 2£ t and %, respectively. 
These are constraint sets defined by linear inequalities that 
involve states and inputs at the same sample time index. 



Motivated by the principles of operator splitting methods 
(see [12] for details and relevant references), we propose 
to split the problem QJ into N + 1 smaller stage-wise 
subproblems that can be solved in parallel. This requires 
breaking the coupling that appears due to the dynamics. We 
introduce a copy of each variable that couples the dynamics 
equations in order to allow such a splitting into subproblems, 
and subsequently impose a consensus constraint on the 
associated complicating variables and their copies. This leads 
to the following equivalent formulation of ([TJ, where the 
complicating variables that are used to perform the splitting 
are clearly highlighted: 



subproblem using 



minimize 



subject to 



-«<' )r R«, ( ' : 



(2a) 



Z f=0 v 

(4 ,wf } )G 9%*%, f = 0,...,iV (2b) 



t = 0,...,N-l 

r (0) - r ■ 
Zt+1 — x t +l 

zt+i=4+i\ t = 0,...,N-l, 



(2c) 

(2d) 
(2e) 
(2f) 



where the subscript t of the decision variables x\ , u\' 
indicates the time index and the superscript (t) denotes the 
group or subproblem where the variable belongs to. Hence, 
each subproblem contains three variables, the current state 
and input as well as a prediction of the state for the next 
time instant. The introduced complicating variable It acts as 
a 'global' variable that brings the local copies x\ f ~ ' and xf' 
in agreement, i.e., z t = x\ ' =x\ '. The time- splitting idea 
is graphically depicted in Figure Q] 
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Fig. 1. The idea of the time splitting algorithm. Whenever the dynamics 

introduce a coupling of variables, they are decoupled via a slack variable 

? -x> -x {t - l) 
z t —x t —x t 



B. The time-splitting algorithm 

In order to use a more compact formulation, we will 
denote the decision variables in (0 corresponding to each 



? _ (J*) ,/) r (0 N 



(3) 



where x t G R 2 "+ m . We also introduce dual variables to deal 
with the consensus equality constraints: 

• w t associated with x\ 1 ' =Zt, t = 1, . . . ,N and 

• v t associated with jq~ ' = zt, t = l,...,N. 

In order to rewrite the finite-time optimal control problem in 
a more compact form, we define the following matrices: 



P t = R t GR (2n+m)x(2n+m)^ (4a) 

_ 

F t = [~A t -B t / ] G R"x(2n+«) 5 ( 4b ) 

G = [ / ] G R«x(2n+m) 5 (4c ) 

Gi = [00/] G R«x( 2 "+^). (4d) 

We use the Alternating Direction Method of Multipliers 
(ADMM) [13], [14] in order to arrive at a solution 
approach that is amenable to parallel implementation. 
The updates involved in the ADMM algorithm include 
forming the augmented Lagrangian of the problem and 
minimizing over the primal variables x t , t = 0,...,N and 
Zt , t = 1, . . . ,iV, followed by updating the dual variables w t 
and Vt, t = l,...,N. The three main steps of the algorithm 
are performed in an iterative fashion and are described 
next in detail. We use k to denote the algorithm's loop 
counter. The termination criterion based on primal and dual 
tolerances are provided and the analytic derivation of the 
formulas in each step is presented in the Appendix. 

Step 1: Solving N+l QP subproblems for the primal 
variables in x t 

Minimization of the augmented Lagrangian over the primal 
variables x t results in Af + 1 stage- wise quadratic programs 
(QPs): 

• For the subproblem associated with the time instant t = 
0, we need to solve 

minimize ( 1 /2)iJP ^o - pv\ T {G\xq - z\) 

+(p/2)||G 1 i -z1l| 2 2 
subject to jco G % 

GoXq = Xi n it 

F x = c 

(5) 
with variable Jco- 

• Similarly, we need to solve the following QPs for all 
the other groups of variables x t , t = 1,...,N— 1: 

minimize ( 1 /2)xJP t x t - pvfjj (G\x t -z k t+l ) 

+(p/2)||G 1 i f -zt + ill! 
-pwf{G x t -~z k t ) 

+ (p/2)\\G x t -z k \\ 2 2 
Xt G &t 
F t x t = c t 



subject to 



(6) 



For the subproblem associated with the final time instant 
t = N, we need to solve 



and we defined the vectors 



minimize 



(l/2)x T N P N x N - pw k J(Gox N - 4) 
+(p/2)||Go%-4||2 

subject to xn G ^v 



with variable jt#. 
The polyhedral sets %, t = 0, . . . ,N are defined as 

^ = ^ x ^ x S% + i C R 2w+m 



(7) 



(8) 



and the variable p > is a parameter of the algorithm. 

Remark 1: Notice that for the time instant N the decision 
variables of the QP actually simplify to xn = xn and ^v = 
3&n, but we keep the same notation for simplicity (and 
without loss of generality). 

Step 2: Averaging 
The update of the 'global' primal variables % , t = 1, . . . ,N is 
derived from a simple quadratic minimization problem, the 
solution of which turns out to be an average of the predicted 
(Xf ') and current (xj ) state 



.**+!. 



.v*+l 






/=!,.. .,#. 



(9) 



This intuitively makes sense, since the global variable can be 
obtained by collecting the local (primal) ones and computing 
the best estimate based on their values. 

Step 3: Dual update 
The dual updates can be expressed as 

w k+l = wf-G if +1 +zf +1 , (10a) 

v* +1 = vf-G 1 if+ 1 +zf + \ f =!,...,#. (10b) 



Termination criterion 

The algorithm terminates when a set of primal and dual 
residuals are bounded by a specified threshold (primal and 
dual tolerances); see [12, §3.2]. The primal and dual residuals 
for the time- splitting algorithm are defined as 



r — A res Xp ri + B res z Dri , 



P n 



■pn ' 



and 



/ = -pALs KS (4-4ri 1 ). 



Pn 



(11) 



(12) 



respectively. 

The termination criterion is activated when 

||r*||2<£ pn , ||/||2<£ dua \ 
where the tolerances e pn and e dual are defined as follows: 

£ P ri = e abs VN2^ 



„dual 



= £ 



(13a) 

+£ rel max{||A res .x pr i||2, ||£ res Vi||2, Ikies W 3b) 
v /(2« + m)(^+l) + e rel ||A^ s v dual || 2 (13c) 
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The residual matrices A res 


and B 


res and the residual vector 


c res are 
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The three update steps described above are fully paral- 
lelizable at each iteration k. Assuming N +1 processors 
are available, then processor Yl t would need to execute the 
following actions for t = 0, . . . ,N — 1: 

1) Receive the estimate z k t+l vf +1 and from neighboring 
processor n f+1 ©. 

2) Compute x k+l . 

3) Receive the estimate S^\ from neighboring processor 
U t i and compute z k+l ©. 

4) Compute w k+l and v k+l 03. 

5) Communicate x k+l to processors n f _i and TL t +i, and 
z k+l , v k+l to processor Yl t \. 

The above scheme suggests that each processor Yl tl t = 
1,...,N — 1 interacts with the two neighboring processors 
n f i and n ?+ i. Processors ITo and 11^ communicate only 
with processors 111 and n#_i, respectively. After updating 
all variables, a gather operation follows in order to compute 
the residuals and check the termination criterion. 

III. THREE-SET SPLITTING QP SOLVER 

A. Motivation 

The time- splitting algorithm presented in the previous sec- 
tion decomposes the centralized finite-time optimal control 
problem so that it can be solved using multiple parallel 
processors. However, the updates for the primal variables 



x t , t = 0, . . . ,iV given in (O, (EJ and ([7} involve solving a 
QP at each iteration of the algorithm. Even though several 
fast interior point solvers exist for this purpose (see e.g., 
[15]), these are mostly suitable for only a limited number of 
variables. Although recently more computationally efficient 
schemes that scale better with the problem size have been 
developed, they are restricted to cases where simple box 
constraints are considered [8], [10], [16]. In order to achieve 
fast computations in an embedded control environment, other 
generic solution methods would be preferred. 

In this section we propose an alternative scheme for 
solving the QPs with general polyhedral constraints that 
arise in the previous section. We propose to perform yet 
another type of splitting approach, which splits the state- 
input variables of the QP in three sets. One set involves 
the variables that appear in the objective function, another 
set includes those that appear in the dynamics equality 
constraints, and the last set contains variables from the 
inequality constraints. In this way, we solve three simpler 
subproblems instead of the single general QP. Since several 
variables are shared among the subproblems, their solutions 
must be in consensus again to ensure consistency. 

An important element of the proposed method is the 
introduction of an extra slack variable, which allows to 
get an analytic solution for the subproblem associated 
with the inequality constraints. Using this variable, the 
projection on any polyhedral set can be practically rewritten 
as a projection onto the nonnegative orthant. Besides this 
feature, the proposed splitting exploits structure in the 
resulting matrices and thus modern numerical linear algebra 
methods can be employed for speeding up the computations. 



B. Problem setup 

We consider a QP of the form 



minimize \x T Mx + q T x + r 
subject to Ax = b 
Hx^h, 



(14) 



with decision variable xGR n , where M G S+, qeR n , r G R, 
A G R mxn , b G R m , H G R pxn , heR p and S n + is the cone of 
positive semidefinite matrices of dimension n. 

In order to apply a variable splitting idea for this problem 
we first replicate all variables appearing in (fl4l three times, 
introducing three different sets for which we must ensure 
consensus. Furthermore, we use a slack variable to remove 
the polyhedral constraint and transform it into a projection 
operation onto the nonnegative orthant. We define the fol- 
lowing sets of variables: 

• First set - objective: x l 

• Second set - equality constraints: x 11 

• Third set - inequality constraints: x m 



• First global variable: z = x 1 = x 11 = x m 

• Second global variable: y = h- Hx^ 

• First set of dual variables: z l ,z ll ,z m associated with 
z = x l = x 11 = x m , respectively. 

• Second set of dual variables: y associated with y = h — 
Hx m 

Using the above sets of variables problem (fT4l can be 
restated in the equivalent form 



minimize 



-x lT Mx l + q T x l + r 



subject to Ax = b 



(15a) 
(15b) 



jn 



y = h-Hx m 1 yhO (15c) 

x l = x ll = x lll = z , (I5d) 

with variables x l ,x J1 ,x m ,z G R n j G R p . The dual variables 
are of dimensions z G R n ,y G R p . 

C. The proposed three-set splitting algorithm 

The proposed algorithm consists of iterative updates to the 
three ADMM steps similarly to the case of the time- splitting 
approach in Section III-B1 namely one for the (local) primal 
variables x l ,x J1 ,x m , one for the (global) primal variables 
z and y and one for the dual variables z l ,z n ,z m and y. 
We provide the algorithm's steps below along with some 
clarifying comments. The analytic derivations are presented 
in the Appendix. 

Step 1: Solving three subproblems for the primal 
variables x 1 ,* 11 ,* 111 

In all three cases we have to solve simple, unconstrained 
QPs. The updates are: 



(xT =(M + piy\p{z k +(?)*) -q) 



pi A T 
A 



(*") 



W k+l 1 



&\k\ 



p(^+(z n ) 

b 



(16) 



(17) 



where v is the dual variable associated with the equality 
constraint Ax 11 = b, and 

(J*)™ = (H T H + l)-\H T (h-y k -f) + 

S + i^f). (18) 

The matrices M + pI and H r // + 7 are symmetric positive 
definite due to the regularization terms. This means that, 
instead of directly inverting the matrices, we can save com- 
putational effort by taking the Cholesky factorization, i.e., 
write the matrix as a product of a lower triangular matrix 
and its transpose (see, e.g., [17]). Furthermore, the matrix 

pi A T 1 
, „ is a KKT matrix and pi >- 0. Hence we can 

A J r 

exploit its structure and use block elimination to solve the 
KKT system (see [18, App. C]). The resulting matrices can 



be pre-factorized and then used in every solve step. The right- 
hand sides are the only parts that change in the update loop. 
Step 2: Averaging and projection 

The update for z is 



1 m 

3 i=l 



Jfc+1 



while for y it is 

y k+l = (h-H(xf 1 ) k+1 -f) . 



(19) 



(20) 



The z-update is an averaging over the three sets of the primal 
variables x 1 ,^ 11 ,^ 111 , while the v-update is the solution of a 
proximal minimization problem (see Appendix), resulting in 
a projection onto the nonnegative orthant, denoted by (•)+. 

Step 3: Dual update 
The update for the dual variables z! is 



a+i 



^k+l 



( ? ) K+, = (?y-W +1 +/ +1 

Similarly, the update for y is 

= y k +y k+l -h+H{x m ) k+l 



i = 1,11,111. (21) 



f+ i 



(22) 



Termination criterion 

The primal and dual tolerances e pn and e dual are given by 

eP ri = e abs ^n~T^^e rel msix{\\A res (x\x l \x m )\\ 2 , 

||5 r esfcv)||2,||c re s||2} 
g dual = £ abs v ^ +£ rel|| A r s( ~I 5 ~II 5 ~III 5 ~ ) || 2 _ 

(23) 
The residual matrices A res and B ms and the residual vector 
c res are 



■^ire.s — 



7 
0/0 
0/ 
0// 



i*rp.s — 



" -/ 


" 


-/ 





-/ 








/ 



and 



c res = (0,0,0,/*), 

where A res e r( 3 "+/>) x3 «, £ res e r( 3 "+/>) ><("+/>) and c res G 
R 3n+P . 

Remark 2: For the QP corresponding to the last sample 
time t = N (|7]>, the algorithm simplifies to splitting into two 
sets (objective and inequality constraints), since there are 
no dynamics equality constraints. The updates and residuals 
follow directly from the more generic case presented above. 

IV. HIERARCHICAL TIME-SPLITTING OPTIMAL 
CONTROL 

It is a natural idea to combine the two splitting algorithms 
(time- splitting and three-set splitting) that were introduced 
in the preceding two sections in order to speed up the 
solution of the finite-time optimal control problem (0Q>. 
This can be accomplished via a nested decomposition 
scheme, where we employ the three-set splitting algorithm 
to solve the QPs Q, © and © appearing in Step 1 of the 
time- splitting algorithm. The idea is graphically depicted in 



Figure 

If we rewrite the generalized inequality constraints appear- 
ing in the problem formulation as 

(4\u^)e% ^H t x t ^h t , t = 0,...,N, (24) 

then the QPs can be written in the form described by (03), 
where we consider the following relations: 

• For xq, Eq. (0: 



M 


= 


P + pG[Gi, 


q 


= 


-pG\ (zi+vi), 


r 


= 


(p/2)z[zi+pv[zi, 


A 


= 


(Go,F ), 


b 


= 


(-*iiiib c o): 


H 


= 


[Ho ] , 


h 


= 


flQ. 



For x t , t = 1,...,JV— 1, Eq. 



M 

q 

r 

A 

b 

H 

h 



= J P f + p(GjG + G[Gi), 

= -p(Gj(z f +w f ) + G[(z m +v m )), 

= (p/2)(zf + izt+i+zjz t )+p(vj +l z t+ i+wjz t ), 

= F t , 

= c t , 

= [H t ] , 

= ht. 



For x N , Eq. (0: 



M 

q 

r 

A 

b 

H 

h 



/V + pGjGo, 

-pGjgv + vPjv), 

o, 
o, 

[H N ] , 
h N . 



For each iteration of the time- splitting algorithm, the three- 
set splitting algorithm runs in an inner loop until it converges. 
The quality of this convergence, i.e., the choice of the primal 
and dual tolerances of the inner loop (|23T > will affect the 
quality of the global solution. 

A method that enables substantial speedup of the algorithm 
is warm starting. Since the three-set splitting algorithm will 
run for every iteration of the time- splitting algorithm, we 
can warm-start each QP with its previous solution. In this 
way, we can achieve a significant reduction in the number 
of iterations needed for the convergence of the inner loop. 



Global Problem 



time splitting 
split 




Global Solution 



Fig. 2. Structural representation of the hierarchical time -splitting solution approach to the finite-time optimal control problem. The outer subproblems at 
t — 0,1,... ,N refer to <|5}, <|6} and 0, respectively. The solutions of the nested subproblems enumerated by I, II and III correspond to ( fTol . JITl and fT8l . 
Notice that only 2 subproblems have to be solved for the last time instant t =N. 



V. NUMERICAL RESULTS 

We consider three, randomly generated, numerical ex- 
amples to illustrate the performance of the algorithm. The 
examples vary in terms of the number of decision vari- 
ables involved. The systems considered are linear and time- 
invariant. We impose constraints on the difference between 
two consecutive states at each time instant of the form 

x t ,i-x t j-i <dx, 

where i/; G R, i = 1, . . . ,n, ,t = 0,...,N and box constraints 
on the inputs, i.e., 

IH|oo < Wmax, t = 0,...,N-l. 

By adjusting the level of the disturbance c t in several time 
instances, we ensure activation of the constraints along the 
horizon. 

For the simulations, we used an Intel Core i7 processor 
running at 1.7 GHz. We compared a C-implementation of 
our algorithm with using CVX [19], a parser-solver that 
uses SDPT3 [20]. For our method, the tolerances for both 
the outer and inner algorithms are set as e pn = 10 4 and 
g duai _ iq 3 j^e parameter p was set after some simple 
tuning. The linear systems appearing in dloV (fTTl and (fT8l 
were solved by first factorizing the matrices off-line, using 



Tim Davis's sparse package [21]-[23] (see also [8]). The 
finite-time optimal control problem was solved only once 
and all the primal and dual variables were initialized at 
zero. However, the inner algorithm was warm- started at 
every iteration of the outer algorithm to the values acquired 
from the previous iteration. No relaxation or any other 
variance of the iterations was used. The numerical results 
are summarized in Table [[] where the computation times are 
reported in ms. 

We can observe that, in the case of the small system, even 
when solving the problem on a single thread, the computation 
times are smaller than those of CVX. As the problem scales, 
the computations have to be parallelized in order to gain 
a significant advantage. More specifically, we expect the 
following speedup factors: 13 and 31 times faster in the case 
of the small problem (for the corresponding tolerances set to 
10 4 and 10 3 respectively). For the medium- sized problem 
the speedups are by a factor of 10.5 and 24.6, and a factor 
of 5.4 and 11.7 for the large-scale problem, respectively. 

In addition, we could observe that the factorization times 
are negligible in all cases, since the matrices being factorized 
are not large. Concerning the three-set splitting algorithm, 
only the average computational times are indicated over all 
iterations required to solve the problem. 





small 


medium 


large 


states n 


10 


20 


50 


inputs m 


10 


10 


40 


horizon length N 


10 


30 


60 


total variables 


220 


900 


5400 


P 


15 


25 


50 


active box constraints 


5 


6 


20 


active inequality constraints 


2 


2 


4 


CVX solve time 


2430 


3529 


19420 


factorization time 


3.18 


9.5 


30 




Tolerance 10~ 4 


3 -set (average) iterations 


21.80 


17 


15.95 


3 -set (average) solve time 


0.75 


1.38 


9.15 


time-split, iterations 


250 


241 


389 


time-split, solve time (single thread) 


1880 


10023 


215888 


time-split, solve time (N threads*) 


188 


334.1 


3598 




Tolerance 10~ 3 


3 -set (average) iterations 


13.14 


13.27 


12.32 


3 -set (average) solve time 


0.49 


1.18 


7.34 


time- split, iterations 


156 


128 


224 


time-split, solve time (single thread) 


780 


4304 


99525 


time-split, solve time (N threads*) 


78 


143.47 


1659 



* estimated parallel computation times 

TABLE I 

Hierarchical time-splitting optimal control: computational 
time results for different size problems. 



VI. CONCLUSIONS 



In this paper, we proposed an algorithm that solves a cen- 
tralized convex finite-time optimal control problem making 
use of operator splitting methods, and, more specifically, 
the Alternating Direction Method of Multipliers. The initial 
problem is split into as many subproblems as the horizon 
length, which then can be solved in parallel. 

The resulting algorithm is composed of three steps, in- 
cluding one where several QPs have to be solved. In this 
respect, we proposed another method, based again on opera- 
tor splitting, that is applicable to QPs of any size, involving 
polyhedral constraints. This algorithm exploits the structure 
of the problem, leading to fast solutions. 

The combination of the proposed algorithms results in a 
nested decomposition scheme for solving the aforementioned 
finite-time optimal control problems over several parallel 
processors. 

Our numerical experiments suggest that the proposed hier- 
archical decomposition approach provides significant speed- 
up in computational time required for medium accuracy 
solutions for the class of problems considered. In our future 
work we intend to perform an even more extensive compar- 
ison with very recent tailor-made computational tools, and 
implement the algorithm on a parallel computing platform to 
obtain more accurate and representative computational time 
measurements. 



APPENDIX 

Derivation of the time -splitting algorithm updates 

We solve the relaxed version of the convex optimization 
problem Q> by formulating the augmented Lagrangian with 
respect to the additional equality constraints, using the ma- 
trices and vectors defined in © and ©. The augmented 
Lagrangian can be written as 

Lp(x t ,x in i t ,w t ,y t ,Zt,Vt):= (25) 

£ I -xfPXt + l % (X t ) J + Vq (x - Xinit) 

+ t t v? +1 (Fx t -c t ) 

t=0 

+ L ( - Pvf+i( G ^t - z t+ i) + -\\Gix t - z t+ i\\l) 
t=o v z ' 

+ t,(-pwJ(GoXt-z t ) + ^\\G x t -Zt\\l). 
t=i K z J 

The dual variables v t e R", t = 0, . . . ,N are associated with 
the equality constraints. We define the indicator function for 
the polyhedral set % c R 2n+m as 

, rn= f o x t e% 

l %\ x t) <y «, otherwise. 

The augmented Lagrangian is minimized (and maximized) 
over the primal (and dual) variables in an iterative manner, 
for each iteration of the algorithm. 

By treating the dynamics and inequality constraints as 
explicit and minimizing with respect to x t , t = 0, . . . ,N, we 
end up with the stage- wise QPs 0, ©, (0. 



The update of It is given by solving 

z k t +l = argmin(-pvf r (G 1 if+ 1 -z f 

7. ^ 



JL " 11/- «*+l s l|2 

+ -\\Gix;z l -z t \\ 2 



zk+\ _ 



pw* T (G x* +1 -z t ) + 9 - ||G if +1 -~z t ||i) 
= argmindll^-Gx^i+ifHi-lllv?!!! 

Go^ + G^-vf-wf 



t = l,...,M26) 



Rearranging the terms in (|251 > and completing the squares 
as before, the dual updates for w t and v t result by solving 

wf+' = argmax(|||w,-Goif +1 +z1 +1 |li-fll*,lli), 

vf +1 = argmaxdlln-G^^+^MH-fllnlll) 



The dual updates are 

■5f+' = vvf-Goif+'+ir 1 , (27) 

v ~*+l = ^-G^+l+^+l, t =\,...,T. (28) 

Substituting _} and <[28> into _} results in wf +1 + 
v* +1 = 0, for all t. Hence I* +1 can be simplified to 

Gojf+' + G.^,' 

This result shows that z f is the variable that enforces consen- 
sus upon Xf and xf~ , the value of which is the average 
of the two. 

Derivation of the three-set splitting algorithm updates 
The augmented Lagrangian for (Q3J can be written as 

L p {x,z,y,zJ,v) := (29) 

ijc^Mjc 1 + q T x l + r + V T (Ax 11 - b) 

in n in 

-p£? T (y-z) + ^£||,'-z||i 



/=! 



2^ 



/=i 
Put 



-pf(h-Hx IU -y) + ^\\h-Hx IU -y\\ 2 2 +I {yeRP+} . 

As before, the dual variable v is associated with the equality 
constraints. We define the indicator function of the nonneg- 
ative orthant for a variable x G R as 

_ _( x x>0 

{xgR + } — + — I oo otherwise. 

The function applies componentwise for vector variables. 

The augmented Lagrangian is a smooth, quadratic 
function with respect to the variables x l ,x J1 ,x m , hence the 
updates can be derived from taking the gradient equal to 
zero. This yields the solutions ([T6T >. (fT7T > and ([T8T >. 



Grouping the terms that involve z from (129b we obtain 
'in 



^ 



J\ k+1 , (j\ k ll2 



= argmin I £||| z - (y)* +i + (?)' 

Ill 



(30) 



Maximizing < f29b over z l , i = I, II, III results in 
(?T = argmax (- f ||^||| + § ||? - (V )*+' + z*+> ||?) . 

Hence, the update for z\i = I, II, III is 

(?)*+' = (g) k - (v) i+1 +/+i, ,- = i,n,m. (3i) 

Similarly, the update for y is 

Taking the sum of the dual variables z\ i = I, II, III from 
([3lT > and substituting into (l30l > simplifies the z-update and we 

obtain 

i m 



3i 



i=\ 



Hence, the global variable z is the average of the three 
estimates of x. 
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